#include "cmatrix.h"

/**
 * CMatrix - creates identity 1x1 matrix
 */
template <typename T>
CMatrix<T>::CMatrix() {
  this->init(1, 1);
}
//---------------------------------------------------------------------------

/**
 * CMatrix - creates identity matrix with given number of rows and columns
 * @param numOfRows number of rows
 * @param numOfCols number of columns
 */
template <typename T>
CMatrix<T>::CMatrix(TIndex numOfRows, TIndex numOfCols) {
  this->init(numOfRows, numOfCols);
}
//---------------------------------------------------------------------------

/**
 * CMatrix - creates matrix from given array of elements
 * @param numOfRows number of rows
 * @param numOfCols number of columns
 * @param items elements array
 */
template <typename T>
CMatrix<T>::CMatrix(TIndex numOfRows, TIndex numOfCols, T items[]) {
  this->init(numOfRows, numOfCols, items);
}
//---------------------------------------------------------------------------

/**
 * CMatrix - creates matrix from vector of CVectors (rows/columns)
 * @param type VECTOR_ROW for rows, VECTOR_COL for columns
 * @param v vector of CVectors
 */
template <typename T>
CMatrix<T>::CMatrix(int type, vector<CVector<T> > v) {
  this->init(type, v);
}
//---------------------------------------------------------------------------

/**
 * CMatrix - creates matrix from array of CVectors (rows/columns)
 * @param type VECTOR_ROW for rows, VECTOR_COL for columns
 * @param numOfVectors number of CVectors in array
 * @param v array of CVectors
 */
template <typename T>
CMatrix<T>::CMatrix(int type, TIndex numOfVectors, CVector<T> v[]) {
  this->init(type, numOfVectors, v);
}
//---------------------------------------------------------------------------

/**
 * CMatrix - copy constructor (deep copy)
 * @param m
 */
template <typename T>
CMatrix<T>::CMatrix(const CMatrix<T> &m) {
  this->copyFrom(m);
}
//---------------------------------------------------------------------------

/**
 * CMatrix - copy constructor (deep copy)
 * @param m
 */
template <typename T>
CMatrix<T>::CMatrix(const CMatrixBase<T> &m) {
  this->copyFrom(m);
}
//---------------------------------------------------------------------------

/**
 * CMatrix - creates matrix by joining matrix and vector
 * @param m matrix
 * @param v vector (will be joined as bottom row or column right column)
 * @param type VECTOR_ROW "v" is row, VECTOR_COL "v" is column
 */
template <typename T>
CMatrix<T>::CMatrix(const CMatrix<T> &m, const CVector<T> &v, int vType) {
  this->init(m, v, vType);
}
//---------------------------------------------------------------------------

/**
 * CMatrix - creates matrix by joining two matrices
 * @param m1 left matrix
 * @param m2 right matrix
 */
template <typename T>
CMatrix<T>::CMatrix(const CMatrix<T> &m1, const CMatrix<T> &m2) {
  this->init(m1, m2);
}
//---------------------------------------------------------------------------

/**
 * ~CMatrix - destructor
 */
template <typename T>
CMatrix<T>::~CMatrix() {
  this->free();
}
//---------------------------------------------------------------------------

/**
 * operator= - assign given matrix to this matrix (does deep copy)
 * @param m
 * @return this
 */
template <typename T>
CMatrix<T>& CMatrix<T>::operator = (const CMatrix<T> &m) {
  this->free();
  this->copyFrom(m);
  return *this;
}
//---------------------------------------------------------------------------

/**
 * operator+= - adds given matrix to this matrix
 * @param m
 * @return this
 */
template <typename T>
CMatrix<T>& CMatrix<T>::operator += (const CMatrix<T> &m) {
  if (numOfRows != m.getNumOfRows() || numOfCols != m.getNumOfCols()) {
    throw CInvalidSizeException(getClassName(), "operator+=");
  }
  return (*this = CMatrixBase<T>::operator+=(m));
}
//---------------------------------------------------------------------------

/**
 * operator-= - subtracts given matrix from this matrix
 * @param m
 * @return this
 */
template <typename T>
CMatrix<T>& CMatrix<T>::operator -= (const CMatrix<T> &m) {
  if (numOfRows != m.getNumOfRows() || numOfCols != m.getNumOfCols()) {
    throw CInvalidSizeException(getClassName(), "operator-=");
  }
  return (*this = CMatrixBase<T>::operator-=(m));
}
//---------------------------------------------------------------------------

/**
 * operator*= - multiplies this matrix by given matrix
 * @param m
 * @return this
 */
template <typename T>
CMatrix<T>& CMatrix<T>::operator *= (const CMatrix<T> &m) {

  // If number of columns of first matrix (factor) differs from number of
  // rows of second matrix (factor), exception
  if (numOfCols != m.getNumOfRows()) {
    throw CInvalidSizeException(getClassName(), "operator*=");
  }
  return (*this = CMatrixBase<T>::operator*=(m));
}
//---------------------------------------------------------------------------

/**
 * operator*= - multiplies this matrix by given number
 * @param num
 * @return this
 */
template <typename T>
CMatrix<T>& CMatrix<T>::operator *= (T num) {
  return (*this = CMatrixBase<T>::operator*=(num));
}
//---------------------------------------------------------------------------

/**
 * operator/= - divides this matrix by given number
 * @param num
 * @return this
 */
template <typename T>
CMatrix<T>& CMatrix<T>::operator /= (T num) {
  return (*this = CMatrixBase<T>::operator/=(num));
}
//---------------------------------------------------------------------------

/**
 * operator+ - adds two matrices
 * @param m
 * @return sum
 */
template <typename T>
CMatrix<T> CMatrix<T>::operator + (const CMatrix<T> &m) const {
  return CMatrix<T>(*this) += m;
}
//---------------------------------------------------------------------------

/**
 * operator- - subtracts two matrices
 * @param m
 * @return difference
 */
template <typename T>
CMatrix<T> CMatrix<T>::operator - (const CMatrix<T> &m) const {
  return CMatrix<T>(*this) -= m;
}
//---------------------------------------------------------------------------

/**
 * operator* - multiply two matrices
 * @param m
 * @return product
 */
template <typename T>
CMatrix<T> CMatrix<T>::operator * (const CMatrix<T> &m) const {
  return CMatrix<T>(*this) *= m;
}
//---------------------------------------------------------------------------

/**
 * operator* - multiply matrix by number
 * @param num
 * @return product matrix
 */
template <typename T>
CMatrix<T> CMatrix<T>::operator * (T num) const {
  return CMatrix<T>(*this) *= num;
}
//---------------------------------------------------------------------------

/**
 * operator/ - divide matrix by number
 * @param num
 * @return divided matrix
 */
template <typename T>
CMatrix<T> CMatrix<T>::operator / (T num) const {
  return CMatrix<T>(*this) /= num;
}
//---------------------------------------------------------------------------

/**
 * setNumOfRows - sets number of rows
 * @param numOfRows number of rows
 */
template <typename T>
void CMatrix<T>::setNumOfRows(TIndex numOfRows) {
  this->numOfRows = numOfRows;
}
//---------------------------------------------------------------------------

/**
 * setNumOfCols - sets number of columns
 * @param numOfCols number of columns
 */
template <typename T>
void CMatrix<T>::setNumOfCols(TIndex numOfCols) {
  this->numOfCols = numOfCols;
}
//---------------------------------------------------------------------------

/**
 * getNumOfRows - returns number of rows
 * @return number of rows
 */
template <typename T>
TIndex CMatrix<T>::getNumOfRows() const {
  return numOfRows;
}
//---------------------------------------------------------------------------

/**
 * getNumOfCols - returns number of columns
 * @return number of columns
 */
template <typename T>
TIndex CMatrix<T>::getNumOfCols() const {
  return numOfCols;
}
//---------------------------------------------------------------------------

/**
 * getRank - calculates rank of matrix (number of LI columns)
 * @return rank of matrix
 */
template <typename T>
TIndex CMatrix<T>::getRank() const {

  // Rank of matrix A (m x n) is h(A) = dim[A*e1, A*e2, ... A*en] where
  // e1, e2, ... en are vectors from standart base En

  // Subtract number of LD columns from total number of columns
  return numOfCols - CMatrix<T>(*this).doGaussElimination();
}
//---------------------------------------------------------------------------

/**
 * isRegular - checks whether is matrix regular
 * @return true if matrix is regular
 */
template <typename T>
bool CMatrix<T>::isRegular() const {

  // Square matrix with size n is regular if n = h(A)

  // If matrix is not square it is not regular, return false
  if (!this->isSquare()) return false;

  // If matrix contains LD columns it is not regular
  return CMatrix<T>(*this).doGaussElimination() > 0 ? false : true;
}
//---------------------------------------------------------------------------

/**
 * getDet - calculates determinant
 * @return determinant
 */
template <typename T>
T CMatrix<T>::getDet() {

  // If matrix is not square, exception
  if (!this->isSquare()) {
    throw CMatrixNotSquareException(getClassName(), "getDet");
  }

  // Temporary variables
  CMatrix<T> backup(*this);
  T det = 1;

  // Iterate through elements on main diagonal and transform matrix to upper
  // triangular
  for (TIndex i = 0; i < numOfRows; i++) {

    // If rows were swapped during pivoting, change sign of determinant
    if (this->doPivoting(i, i)) det *= -1;
    
    // Multiply determinant with actual diagonal element, if element is zero
    // => determinant is zero, break cycle
    if ((det *= this->items[i][i]) == 0) break;
    
    // Eliminate column under actual diagonal element, if any row was
    // multiplied/divided by some number then determinant must be
    // divided/multiplied by same number
    det *= this->elimColUnder(i, i);
  }

  // Matrix is upper triangular all eigen values lies on main diagonal =>
  // determinant is equal to product of diagonal elements multiplied by "det"
  // because all diagonal elements are one => determinant is "det" * 1
  *this = backup;
  return det;
}
//---------------------------------------------------------------------------

/**
 * invert - inverts matrix (calculates inverse matrix)
 */
template <typename T>
void CMatrix<T>::invert() {

  // If matrix is not square, exception
  if (!this->isSquare()) {
    throw CMatrixNotSquareException(getClassName(), "invert");
  }

  // Create temporary matrix (join this matrix from right by identity matrix)
  CMatrix<T> tmpMatrix = CMatrix<T>(*this, CMatrix<T>(numOfRows, numOfCols));

  // Perform full elimination of temporary matrix, if there are some LD
  // columns => original matrix wasn't regular => cannot be inverted,
  // exception
  if (!tmpMatrix.doGaussJordanElimination()) {
    throw CMatrixNotRegException(getClassName(), "invert");
  }

  // Get resulting inverse matrix from temporary matrix and copy it to this
  *this = tmpMatrix.subMatrix(0, numOfCols, numOfRows, numOfCols * 2);
}
//---------------------------------------------------------------------------

/**
 * getBase - gets base of vector subspace, made from LI columns of matrix
 * @return submatrix
 */
template <typename T>
CMatrix<T> CMatrix<T>::getBase() {

  // Temporary matrix (copy of this) and number of LI columns
  CMatrix<T> tmp(*this);
  TIndex tmpNumOfCols = numOfCols - tmp.doGaussElimination();

  // If number of LI columns is equal to total number of columns, base is
  // this
  if (tmpNumOfCols == numOfCols) return *this;

  // Resulting LI columns
  CVector<T> tmpCols[tmpNumOfCols];

  // Go through diagonal elements of temporaty matrix (in upper triangular
  // form), calculate indices of main (LI) columns and get these columns from
  // this
  TIndex i = 0, j = 0;
  while (i < numOfRows && j < numOfCols) {
    if (!tmp.getItemsPtr()[i][j]) {
      j++;
      continue;
    }
    tmpCols[i++] = this->getCol(j++);
  }
  return CMatrix<T>(this->VECTOR_COL, tmpNumOfCols, tmpCols);
}
//---------------------------------------------------------------------------

/**
 * subMatrix - gets submatrix from matrix
 * @param sr starting row index
 * @param sc starting column index
 * @param er ending row index
 * @param ec ending column index
 * @return submatrix
 */
template <typename T>
CMatrix<T> CMatrix<T>::subMatrix(TIndex sr, TIndex sc, TIndex er, TIndex ec) {

  // Check whether indices are valid (pointing to matrix)
  if (sr >= numOfRows || er > numOfRows || sc >= numOfCols || ec > numOfCols) {
    throw CIndexOutOfBoundsException(getClassName(), "subMatrix");
  }

  // Number of rows, columns and elements of submatrix
  TIndex tmpNumOfRows = er - sr;
  TIndex tmpNumOfCols = ec - sc;
  T tmpItems[tmpNumOfRows * tmpNumOfCols];

  // Create elements of submatrix
  for (TIndex i = sr; i < er; i++) {
    for (TIndex j = sc; j < ec; j++) {
      tmpItems[i * (tmpNumOfCols - sr) + (j - sc)] = this->items[i][j];
    }
  }
  return CMatrix<T>(tmpNumOfRows, tmpNumOfCols, tmpItems);
}
//---------------------------------------------------------------------------

/**
 * mul - multiplies matrix by vector
 * @param m matrix
 * @param v vector
 * @return product
 */
template <typename T>
CVector<T> CMatrix<T>::mul(const CMatrix<T> &m, const CVector<T> &v) {

  // If number of matrix columns differs from vector size, exception
  if (m.getNumOfCols() != v.getSize()) {
    throw CInvalidSizeException(getClassName(), "mul");
  }
  return CMatrixBase<T>::mul(m, v);
}
//---------------------------------------------------------------------------

/**
 * mul - multiplies vector by matrix
 * @param v vector
 * @param m matrix
 * @return product
 */
template <typename T>
CVector<T> CMatrix<T>::mul(const CVector<T> &v, const CMatrix<T> &m) {

  // If number of matrix rows differs from vector size, exception
  if (v.getSize() != m.getNumOfRows() ) {
    throw CInvalidSizeException(getClassName(), "mul");
  }
  return CMatrixBase<T>::mul(v, m);
}
//---------------------------------------------------------------------------

template class CMatrix<char>;
template class CMatrix<short>;
template class CMatrix<int>;
template class CMatrix<long>;
template class CMatrix<float>;
template class CMatrix<double>;
template class CMatrix<CBigInt>;
template class CMatrix<CFraction>;
